Morphological divergence in fragmented populations: No evidence of fitness benefit

The code below was used to produce all figures and tables used in the paper.

Data files ADD_I and ADD_R are available on the github page.

Study Design

R setup

Load all packages

# libraries
library(readr)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(gganimate)
library(gifski)
library(dplyr)
library(tidyr)
library(ggthemes)
library(ggpubr)
library(lme4)
library(ggeffects)
library(lmerTest)
library(car)
library(ggplot2)
library(survival)
library(coxme)
library(survminer)
library(ggfortify)
library(broom)
library(ehahelper)
library(dena)
library(tidyr)

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/London
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] dena_0.1.0         ehahelper_0.3.9999 broom_1.0.4        ggfortify_0.4.16  
 [5] survminer_0.4.9    coxme_2.2-18.1     bdsmatrix_1.3-6    survival_3.5-5    
 [9] car_3.1-2          carData_3.0-5      lmerTest_3.1-3     ggeffects_1.2.2   
[13] lme4_1.1-33        Matrix_1.5-4       ggpubr_0.6.0       ggthemes_4.2.4    
[17] tidyr_1.3.0        gifski_1.6.6-1     gganimate_1.0.8    ggplot2_3.4.2     
[21] sjlabelled_1.2.0   sjmisc_2.8.9       sjPlot_2.8.14      readr_2.1.4       
[25] data.table_1.14.8  dplyr_1.1.2       

loaded via a namespace (and not attached):
 [1] gridExtra_2.3       rlang_1.1.1         magrittr_2.0.3      compiler_4.3.0     
 [5] mgcv_1.8-42         vctrs_0.6.2         stringr_1.5.0       pkgconfig_2.0.3    
 [9] crayon_1.5.2        fastmap_1.1.1       backports_1.4.1     labeling_0.4.2     
[13] KMsurv_0.1-5        effectsize_0.8.3    utf8_1.2.3          rmarkdown_2.21     
[17] tzdb_0.3.0          haven_2.5.2         nloptr_2.0.3        purrr_1.0.1        
[21] bit_4.0.5           xfun_0.39           cachem_1.0.8        jsonlite_1.8.4     
[25] progress_1.2.2      tweenr_2.0.2        parallel_4.3.0      prettyunits_1.1.1  
[29] R6_2.5.1            bslib_0.4.2         stringi_1.7.12      pkgload_1.3.2      
[33] boot_1.3-28.1       jquerylib_0.1.4     numDeriv_2016.8-1.1 estimability_1.4.1 
[37] Rcpp_1.0.10         knitr_1.42          modelr_0.1.11       zoo_1.8-12         
[41] parameters_0.21.0   splines_4.3.0       tidyselect_1.2.0    yaml_2.3.7         
[45] rstudioapi_0.14     abind_1.4-5         lattice_0.21-8      tibble_3.2.1       
[49] withr_2.5.0         bayestestR_0.13.1   evaluate_0.20       survMisc_0.5.6     
[53] pillar_1.9.0        insight_0.19.1      plotly_4.10.1       generics_0.1.3     
[57] vroom_1.6.3         hms_1.1.3           munsell_0.5.0       scales_1.2.1       
[61] minqa_1.2.5         xtable_1.8-4        glue_1.6.2          emmeans_1.8.5      
[65] lazyeval_0.2.2      tools_4.3.0         ggsignif_0.6.4      forcats_1.0.0      
[69] mvtnorm_1.1-3       cowplot_1.1.1       grid_4.3.0          datawizard_0.7.1   
[73] colorspace_2.1-0    nlme_3.1-162        eha_2.10.3          performance_0.10.3 
[77] cli_3.6.1           km.ci_0.5-6         fansi_1.0.4         viridisLite_0.4.2  
[81] sjstats_0.18.2      gtable_0.3.3        rstatix_0.7.2       sass_0.4.6         
[85] digest_0.6.31       htmlwidgets_1.6.2   farver_2.1.1        htmltools_0.5.5    
[89] lifecycle_1.0.3     httr_1.4.5          bit64_4.0.5         MASS_7.3-60        

Morphometric and Island Comparison

Preparing Data

# read the data file "ADD_I"
ADD_I <- read_csv("Tables/ADD_I.csv", col_types = cols(Island = col_factor(levels = c("CN","AR", "CE", "DS", "FR"))))
# create the earliest born individuals on each island
Isfy <- ADD_I %>% select(founderyear,Island) %>% filter(!duplicated(Island))
# create a data frame with means and sd of each morphometric
ADD_mean <- ADD_I %>% group_by(Island,BirthYear) %>% summarise(RTsd = sd(RightTarsus,na.rm=TRUE),WLsd = sd(WingLength,na.rm=TRUE),BMsd = sd(BodyMass,na.rm=TRUE),HBsd=sd(HeadBill,na.rm=TRUE),RightTarsus=mean(RightTarsus,na.rm=TRUE),WingLength=mean(WingLength,na.rm=TRUE),BodyMass=mean(BodyMass,na.rm=TRUE),HeadBill=mean(HeadBill,na.rm=TRUE))

`summarise()` has grouped output by 'Island'. You can override using the `.groups` argument.

Linear Mixed Effect Models

Tarsus Length

####RightTarsus ----
IRTlmer <- lmer(RightTarsus ~  Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IRTlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RightTarsus ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 10511

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.0560 -0.4445  0.0052  0.4324  5.5726 

Random effects:
 Groups   Name        Variance  Std.Dev.
 BirdID   (Intercept) 0.3048654 0.55215 
 Observer (Intercept) 0.0562847 0.23724 
 Ageclass (Intercept) 0.0009875 0.03142 
 Residual             0.1407944 0.37523 
Number of obs: 6154, groups:  BirdID, 3258; Observer, 82; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     2.417e+01  6.083e-02  2.713e+02 397.385  < 2e-16 ***
IslandAR       -7.215e-02  1.686e-01  3.839e+03  -0.428 0.668783    
IslandCE       -2.510e-01  9.829e-02  3.925e+03  -2.553 0.010704 *  
IslandDS        2.519e-01  8.080e-02  3.408e+03   3.117 0.001842 ** 
IslandFR        1.869e-01  1.236e-01  2.500e+03   1.512 0.130545    
fYear           8.838e-03  2.530e-03  2.174e+03   3.493 0.000488 ***
SexMale         1.596e+00  2.244e-02  3.211e+03  71.130  < 2e-16 ***
IslandAR:fYear -6.760e-04  1.024e-02  3.539e+03  -0.066 0.947388    
IslandCE:fYear  1.764e-03  6.196e-03  4.079e+03   0.285 0.775832    
IslandDS:fYear -2.910e-03  5.840e-03  4.085e+03  -0.498 0.618312    
IslandFR:fYear  1.554e-03  1.756e-02  1.981e+03   0.089 0.929486    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.163                                                               
IslandCE    -0.309  0.111                                                        
IslandDS    -0.473  0.113  0.218                                                 
IslandFR    -0.331  0.081  0.150  0.225                                          
fYear       -0.779  0.188  0.339  0.504  0.366                                   
SexMale     -0.177 -0.029  0.008 -0.033 -0.024 -0.016                            
IslndAR:fYr  0.149 -0.922 -0.096 -0.101 -0.076 -0.200  0.026                     
IslndCE:fYr  0.261 -0.094 -0.915 -0.173 -0.126 -0.323 -0.010  0.097              
IslndDS:fYr  0.314 -0.065 -0.138 -0.849 -0.208 -0.379  0.008  0.069  0.123       
IslndFR:fYr  0.102 -0.025 -0.050 -0.108 -0.850 -0.143  0.020  0.029  0.048  0.215

#predicting the model 
IRTlmerdata <- ggpredict(IRTlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,RightTarsus=predicted,Island=group,Sex=facet)
IRTlmerdata2 <- merge(IRTlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(RightTarsus = mean(RightTarsus, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model (shown as ggarrange below)
IRTmod <- ggplot(IRTlmerdata2, aes(x = BirthYear, y = RightTarsus, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=RightTarsus,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=RightTarsus-RTsd,ymax=RightTarsus+RTsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "loess", se = FALSE) + 
  geom_ribbon(aes(ymin=RightTarsus-std.error,ymax=RightTarsus+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Tarsus Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Body Mass

#### Body Mass ----
IBMlmer <- lmer(BodyMass ~  Island*fYear + Sex  + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IBMlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BodyMass ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 16677.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6968 -0.5705 -0.0271  0.5440  4.2625 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.24508  0.4951  
 Observer (Intercept) 0.05818  0.2412  
 Ageclass (Intercept) 0.14657  0.3828  
 Residual             0.53384  0.7306  
Number of obs: 6554, groups:  BirdID, 3398; Observer, 91; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     1.438e+01  1.846e-01  5.339e+00  77.909 2.29e-09 ***
IslandAR        1.189e+00  1.570e-01  3.120e+03   7.577 4.63e-14 ***
IslandCE        2.792e-01  1.214e-01  3.666e+03   2.299  0.02158 *  
IslandDS        7.385e-01  1.037e-01  1.096e+03   7.122 1.92e-12 ***
IslandFR        5.264e-01  1.668e-01  1.471e+03   3.156  0.00163 ** 
fYear          -2.623e-03  3.107e-03  5.626e+02  -0.844  0.39893    
SexMale         1.433e+00  2.635e-02  2.907e+03  54.405  < 2e-16 ***
IslandAR:fYear -7.870e-02  1.043e-02  2.346e+03  -7.546 6.38e-14 ***
IslandCE:fYear -3.708e-02  7.935e-03  4.479e+03  -4.673 3.05e-06 ***
IslandDS:fYear -5.069e-02  7.678e-03  2.097e+03  -6.601 5.14e-11 ***
IslandFR:fYear -5.400e-02  2.379e-02  1.097e+03  -2.270  0.02343 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.071                                                               
IslandCE    -0.107  0.098                                                        
IslandDS    -0.176  0.107  0.184                                                 
IslandFR    -0.113  0.071  0.120  0.187                                          
fYear       -0.307  0.188  0.294  0.462  0.320                                   
SexMale     -0.068 -0.016  0.007 -0.035 -0.027 -0.017                            
IslndAR:fYr  0.065 -0.869 -0.082 -0.095 -0.067 -0.205  0.016                     
IslndCE:fYr  0.092 -0.086 -0.910 -0.144 -0.097 -0.283 -0.007  0.092              
IslndDS:fYr  0.117 -0.064 -0.114 -0.847 -0.181 -0.345  0.015  0.069  0.098       
IslndFR:fYr  0.035 -0.024 -0.039 -0.097 -0.868 -0.127  0.023  0.027  0.035  0.201

#predicting the model
IBMlmerdata <- ggpredict(IBMlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,BodyMass=predicted,Island=group,Sex=facet)
IBMlmerdata2 <- merge(IBMlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)  %>% group_by(Island,BirthYear) %>% summarise(BodyMass = mean(BodyMass, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IBMmod <- ggplot(IBMlmerdata2, aes(x = BirthYear, y = BodyMass, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=BodyMass,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=BodyMass-BMsd,ymax=BodyMass+BMsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=BodyMass-std.error,ymax=BodyMass+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Body Mass (g)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Wing Length

#### Wing Length ----
IWLlmer <- lmer(WingLength ~  Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)

Warning: Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)

summary(IWLlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: WingLength ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 23050.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8213 -0.5350  0.0161  0.5335  4.0139 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.95940  0.9795  
 Observer (Intercept) 0.40157  0.6337  
 Ageclass (Intercept) 0.08901  0.2983  
 Residual             1.26489  1.1247  
Number of obs: 6510, groups:  BirdID, 3408; Observer, 91; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     6.592e+01  1.905e-01  1.616e+01 346.012  < 2e-16 ***
IslandAR        2.490e-01  2.795e-01  3.717e+03   0.891  0.37296    
IslandCE       -1.869e-01  2.098e-01  4.116e+03  -0.891  0.37310    
IslandDS       -5.162e-01  1.865e-01  2.916e+03  -2.767  0.00569 ** 
IslandFR        1.871e-01  2.876e-01  2.198e+03   0.650  0.51549    
fYear          -2.900e-02  5.796e-03  1.679e+03  -5.004 6.22e-07 ***
SexMale         2.589e+00  4.624e-02  3.176e+03  55.991  < 2e-16 ***
IslandAR:fYear -4.495e-02  1.871e-02  3.681e+03  -2.403  0.01632 *  
IslandCE:fYear  2.268e-03  1.354e-02  4.598e+03   0.168  0.86696    
IslandDS:fYear  4.375e-02  1.350e-02  4.156e+03   3.241  0.00120 ** 
IslandFR:fYear  2.893e-02  4.134e-02  2.172e+03   0.700  0.48407    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.123                                                               
IslandCE    -0.198  0.099                                                        
IslandDS    -0.315  0.101  0.190                                                 
IslandFR    -0.225  0.072  0.131  0.197                                          
fYear       -0.551  0.177  0.302  0.464  0.344                                   
SexMale     -0.118 -0.015  0.005 -0.033 -0.026 -0.015                            
IslndAR:fYr  0.115 -0.873 -0.083 -0.092 -0.070 -0.200  0.013                     
IslndCE:fYr  0.174 -0.088 -0.910 -0.152 -0.112 -0.299 -0.007  0.096              
IslndDS:fYr  0.218 -0.061 -0.121 -0.851 -0.195 -0.359  0.012  0.068  0.107       
IslndFR:fYr  0.074 -0.023 -0.043 -0.099 -0.857 -0.135  0.024  0.028  0.040  0.217
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)

#predicting the model
IWLlmerdata <- ggpredict(IWLlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,WingLength=predicted,Island=group,Sex=facet)
IWLlmerdata2 <- merge(IWLlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11)) %>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear)%>% summarise(WingLength = mean(WingLength, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IWLmod <- ggplot(IWLlmerdata2, aes(x = BirthYear, y = WingLength, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=WingLength,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=WingLength-WLsd,ymax=WingLength+WLsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=WingLength-std.error,ymax=WingLength+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Wing Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Head + Bill Length

#### Head + Bill ----
IHBlmer <- lmer(HeadBill ~  Island*fYear + Sex  + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IHBlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: HeadBill ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 9601.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1841 -0.4234  0.0279  0.4642  5.4733 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.22985  0.4794  
 Observer (Intercept) 0.04913  0.2217  
 Ageclass (Intercept) 0.04152  0.2038  
 Residual             0.14507  0.3809  
Number of obs: 5937, groups:  BirdID, 3120; Observer, 80; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     3.588e+01  1.085e-01  7.815e+00 330.737   <2e-16 ***
IslandAR        4.153e-01  2.063e-01  1.204e+03   2.013   0.0443 *  
IslandCE       -1.041e-01  9.292e-02  3.741e+03  -1.121   0.2624    
IslandDS        8.141e-02  7.652e-02  3.252e+03   1.064   0.2874    
IslandFR       -7.000e-02  1.160e-01  2.381e+03  -0.603   0.5464    
fYear          -1.237e-03  2.492e-03  2.573e+03  -0.496   0.6196    
SexMale         8.420e-01  2.075e-02  3.016e+03  40.579   <2e-16 ***
IslandAR:fYear -2.977e-02  1.559e-02  5.975e+02  -1.910   0.0566 .  
IslandCE:fYear -8.098e-03  5.841e-03  3.978e+03  -1.386   0.1657    
IslandDS:fYear  1.404e-03  5.456e-03  3.918e+03   0.257   0.7969    
IslandFR:fYear  1.196e-02  1.638e-02  1.925e+03   0.730   0.4655    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.066                                                               
IslandCE    -0.177  0.079                                                        
IslandDS    -0.277  0.088  0.233                                                 
IslandFR    -0.197  0.067  0.163  0.243                                          
fYear       -0.444  0.155  0.351  0.525  0.386                                   
SexMale     -0.089 -0.017 -0.001 -0.036 -0.029 -0.021                            
IslndAR:fYr  0.046 -0.944 -0.049 -0.062 -0.051 -0.133  0.009                     
IslndCE:fYr  0.152 -0.064 -0.918 -0.188 -0.138 -0.333 -0.002  0.047              
IslndDS:fYr  0.187 -0.049 -0.147 -0.851 -0.220 -0.393  0.012  0.041  0.131       
IslndFR:fYr  0.062 -0.021 -0.054 -0.113 -0.848 -0.148  0.023  0.019  0.050  0.218

#predicting the model
IHBlmerdata <- ggpredict(IHBlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,HeadBill=predicted,Island=group,SexEstimate=facet)
IHBlmerdata2 <- merge(IHBlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 22)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)%>% group_by(Island,BirthYear)%>% summarise(HeadBill = mean(HeadBill, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IHBmod <- ggplot(IHBlmerdata2, aes(x = BirthYear, y = HeadBill, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=HeadBill,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=HeadBill-HBsd,ymax=HeadBill+HBsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=HeadBill-std.error,ymax=HeadBill+std.error),alpha=0.25) +
  labs(x = "Birth Year", y= "Head + Bill Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Arrange and plot

Imods <- ggarrange(IRTmod,IBMmod,IWLmod,IHBmod, common.legend=TRUE, labels = c("A","B","C","D"))

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'

Imods

Lifespan and Reproduction

Preparing data

ADD_R <- read_csv("Tables/ADD_R.csv", col_types = cols(SexEstimate = col_factor()))
ADD_Real <- as.data.frame(ADD_R)
ADD_Real <- ADD_Real %>% select(RightTarsus,BodyMass,WingLength,HeadBill,BirdID,SexEstimate,BirthYear,Lifespan,SurvStatus,ReproductiveOutput) 
ADD_tidy <- ADD_Real %>% filter(!is.na(RightTarsus) & !is.na(BodyMass)  & !is.na(WingLength)  & !is.na(HeadBill) & !is.na(SexEstimate)  & !is.na(Lifespan) ) %>% mutate(BodyCondition = BodyMass/RightTarsus) %>% mutate(SexEstimate = as.factor(SexEstimate))

Survival analysis with coxme

Model comparisons

###Lifespan analysis with coxme----

lscox0 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ BodyMass + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox1 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox2 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus*poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox3 <- coxme(Surv(Lifespan,SurvStatus) ~ poly(BodyCondition,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox4 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2)*SexEstimate + WingLength + HeadBill + (1|BirthYear), data= ADD_tidy)

AIC(lscox0,lscox1,lscox2,lscox3,lscox4) #quadratic body mass is best

BIC(lscox0,lscox1,lscox2,lscox3,lscox4)

# lscox4 was selected based on the best AIC and BIC values

Predicting and plotting coxme

####predict lifespan----
ADD_round <- ADD_R %>%mutate(across(c(RightTarsus, BodyMass, HeadBill, WingLength), round, 0)) %>% distinct(RightTarsus, BodyMass, HeadBill, WingLength,SexEstimate) %>% filter(!is.na(RightTarsus) & !is.na(BodyMass)  & !is.na(WingLength)  & !is.na(HeadBill) & !is.na(SexEstimate))

ADD_end <- bind_rows(ADD_tidy,ADD_round)

rs <- predict_coxme(lscox4,newdata = ADD_end,type="risk", se.fit=TRUE)
rsend <- ADD_end %>% add_columns(rr = rs$fit[,1], se.fit = rs$se.fit[,1])
ADD_tail <- rsend %>% tail(nrow(ADD_round)) %>% mutate(SexEstimate = as.factor(SexEstimate)) %>% mutate(low = exp(log(rr) - 1.96 * se.fit),high = exp(log(rr) + 1.96 * se.fit)) %>% group_by(BodyMass,SexEstimate) %>% summarise(mean_rr=mean(rr),mean_se.fit=mean(se.fit), mean_low=mean(low),mean_high=mean(high))

`summarise()` has grouped output by 'BodyMass'. You can override using the `.groups` argument.

coxmeplot <- ggplot(ADD_tail,aes(x=BodyMass,y=mean_rr)) + 
  geom_line(aes(colour=SexEstimate)) +
  geom_ribbon(aes(fill=SexEstimate,ymin=mean_rr-mean_se.fit,ymax=mean_rr+mean_se.fit), alpha = 0.25) + 
  labs(x = "Body Mass (g)", y = "Mortality Risk Score") +
  guides(fill=guide_legend(title="Sex"), colour=guide_legend(title="Sex")) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind(labels = c("Female", "Male")) +
  scale_fill_colorblind(labels = c("Female", "Male"))

Reproductive Output using GLM

Model comparisons

ROglmF <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill  ,data = ADD_tidy, family = "poisson")
ROglmR <- glm(ReproductiveOutput ~ poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm0 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill  + Lifespan + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm1 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength  + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm2 <- glm(ReproductiveOutput ~ poly(RightTarsus,2)+ BodyMass + WingLength + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm3 <- glm(ReproductiveOutput ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm4 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + poly(WingLength,2) + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm5 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + poly(HeadBill,2)  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")

AIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # 4 best, 1 second

BIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # 1 best 

vif(ROglm1)

                      GVIF Df GVIF^(1/(2*Df))
RightTarsus       3.252977  1        1.803601
BodyMass          2.553710  1        1.598033
WingLength        2.550319  1        1.596972
HeadBill          2.274822  1        1.508251
poly(Lifespan, 2) 1.212419  2        1.049333
SexEstimate       3.638709  1        1.907540
BirthYear         1.266770  1        1.125509

summary(ROglm1)

Call:
glm(formula = ReproductiveOutput ~ RightTarsus + BodyMass + WingLength + 
    HeadBill + poly(Lifespan, 2) + SexEstimate + BirthYear, family = "poisson", 
    data = ADD_tidy)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        55.943199   6.971746   8.024 1.02e-15 ***
RightTarsus        -0.004860   0.032759  -0.148  0.88207    
BodyMass            0.075016   0.028836   2.602  0.00928 ** 
WingLength          0.015426   0.015372   1.004  0.31561    
HeadBill           -0.002388   0.038924  -0.061  0.95108    
poly(Lifespan, 2)1 23.108050   1.113788  20.747  < 2e-16 ***
poly(Lifespan, 2)2 -5.454150   0.702063  -7.769 7.93e-15 ***
SexEstimate1       -0.042798   0.066926  -0.639  0.52251    
BirthYear          -0.028461   0.003275  -8.692  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2018.43  on 736  degrees of freedom
Residual deviance:  879.58  on 728  degrees of freedom
  (629 observations deleted due to missingness)
AIC: 3146.9

Number of Fisher Scoring iterations: 5

Predicting and plotting GLM

ROplotdata <- ADD_tidy %>% mutate(BodyMass = case_when(BodyMass > 16.5 ~ "16.5", BodyMass < 14.5 ~ "14.5",BodyMass < 16.5 & BodyMass > 14.5 ~ "15.5")) %>% filter(!is.na(BodyMass))
ROglmdat <- ggpredict(ROglm1, terms = c("Lifespan","BodyMass")) %>% rename(Lifespan=x,ReproductiveOutput=predicted,BodyMass=group)
ROglmplot <- ggplot(ROglmdat, aes(x=Lifespan,y=ReproductiveOutput)) + 
  geom_point(data = ROplotdata, mapping=aes(x=Lifespan,y=ReproductiveOutput,colour=BodyMass),position="jitter") +
  geom_smooth(aes(colour=BodyMass), se = FALSE) + 
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high,fill=BodyMass), alpha = 0.25) +
  labs(x = "Lifespan (years)", y= "Total Number of Offspring") +
  theme_tufte(base_size = 15, base_family = "Arial") +
  scale_color_colorblind() + scale_fill_colorblind()

Arrange and plot

fitplot <- ggarrange(coxmeplot,ROglmplot,  labels = c("A","B"))

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

fitplot

Thank you

Thank you for going through the code and supporting open research.